Thermalization in a coherently driven ensemble of two-level systems 
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We study the coherent quantum evolution of a closed and driven mesoscopic chain of two-level 
systems that interact via the van-der-Waals interaction in their excited state. The Hamiltonian 
consists of a part corresponding to a classical lattice gas and an off-diagonal driving term without 
classical counterpart. We show that in a certain parameter range the latter leads to a thermalization 
of the system with respect to observables of the classical lattice gas such as the interaction energy 
and particle number distribution. We investigate the evolution of the system into this thermal 
state and discuss how to determine the corresponding temperature. Our findings can be applied to 
understand thermalization in strongly interacting systems of laser-driven Rydberg atoms, ions or 
polar molecules. 
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One of the most intriguing questions in many-body 
physics is how and under what conditions a closed quan- 
tum system evolves into thermal equilibrium. This prob- 
lem has recently received much attention and the ap- 
proach of a closed system to equilibrium (or the absence 
of it) has been subject of a number of experimental and 
theoretical efforts [l|-|8[. The former have been predomi- 
nantly driven by the progress in ultracold atomic physics 
which allow to study coherent many-body dynamics in a 
very clean and decoherence-free environment Q. 

In this work we investigate the thermalization of a 
mesoscopic ensemble of strongly and long-range interact- 
ing (laser-) driven two-level systems represented by atoms 
@i 3 Hi; ions [l^l or polar molecules fisl] confined in 
a deep lattice. The external dynamics is frozen out and 
the Hamiltonian of the internal dynamics can be decom- 
posed into a classical part -//lg corresponding to a lattice 
gas and a quantum part Hq originating from the driv- 
ing. We show that in a certain parameter regime the 
system reaches a steady state after an initial transient 
period. Here observables of the classical lattice gas, such 
as the mean number of excited particles and the inter- 
action energy, have thermalized according to a canonical 
distribution with of temperature T. We show that T 
is proportional to the strength of the laser driving, and 
outline two ways for its experimental determination. The 
first method relies on measuring the statistics of the inter- 
action energies. The second is based on determining the 
response of the number of excited particles to a change 
of the chemical potential. Both methods are applicable 
to recent experiments with ultracold and highly excited 
atoms 
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We aim for our analysis to stay close to experimen- 
tally realizable conditions. We therefore consider initial 
conditions which can be achieved and observables which 
can be accessed in actual experiments. We hope that our 
results here will motivate new experimental research to 
shed light onto the general issue of thermalization. 

Our system consists of a periodic chain of L driven two- 



level systems which - once excited - interact via the van- 
der-Waals interaction which decays with the inverse sixth 
power of the interparticle distance. The Hamiltonian is 
given hy H = Hq + Hi^q - fi Ylk=i "-fc "^it^i 



He 



fe=i 



L 

E 



m 



(1) 



The ground state \g)f, of each particle is coupled to its 



excited state |e),, via the operator a. 



(k) 



|e),(5l + l5),(e| 



at a rate (Rabi frequency) This form is typical for a 
laser-atom interaction treated within the rotating wave 
approximation. The operator = |e);;(e| projects on 
the excited state of the fc-th particle. The energy scale 
of the interaction is given ^ayV = C^/a^ where Cq is the 
van-der-Waals coefficient [l3, and a is the (regular) 
spacing between the two-level systems. The quantity 
corresponds to the chemical potential which in practice 
can be adjusted by varying the detuning of the laser that 
drives the transition: \g) — |e). 

Throughout this paper we focus on the coherent evo- 
lution of the system from the initial (product) state 
|init) = ®k=i Iff)*;' ^^^^ state is the natural start- 
ing point from an experimental perspective. We are in- 
terested in the approach of the system to equilibrium. 
Since the system is closed, and evolves coherently, equi- 
libration is not expected to occur with respect to any 
arbitrary observable Here we focus on observables 
of the lattice gas represented through /Jlg • Such 'clas- 
sical' observables are the total number of excited par- 
ticles, N = rife, or the density-density correlations, 

gnra = L"^ {Umnn) / {Nf. 

In Fig. [Ua) we show the time-evolution of these lat- 
tice gas quantities for the parameter choice V = bfl for 
which equilibration seems to occur. The observed be- 
havior is the same for all lattice gas observables shown: 
The initial large amplitude oscillations diminish in the 
interval 10^ < fit < 10^, and for even longer times they 
are drastically reduced — the subsystem probed by these 
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FIG. 1: (a) Temporal evolution of the mean excitation num- 
ber (A'^) (t) and the correlation functions (?i2(i) and giait) 
{L = 20, V = 5fl, = 0). All quantities reach their steady- 
state - characterized by small amplitude oscillations - in the 
interval 10^ < fit < 10^ . (b) Expectation value of the number 
of excited particles (A'') (blue) and its fluctuations AA'^ (red), 
studied over the time interval 47 < Sit < 50, as a function of 
the interaction strength. The thin lines correspond to 20 in- 
stantaneous snapshots taken during that period. Thick solid 
lines represent the average over that time interval. For very 
small V no relaxation into a steady state is occurring as large 
oscillations of the (A'^) and AA'" are present. For large values 
of V the instantaneous values of (A'^) and AA^ differ only very 
little from the averaged ones. Moreover, the fluctuations AA'' 
are small. We refer to this as the thermal regime. 



observables appears to have reached a steady state. Sim- 
ilar behavior has been reported for other systems such as 
hard core bosons confined to a lattice [l^ . In the follow- 
ing we will mainly focus on the expectation value of N 
with respect to which the system reaches a steady state 
a,t fit Ri 50 for the parameters of Fig. [TJa). All times 
considered in this work are smaller than the quantum 
mechanical revival time. 

Our aim is now to analyze why and under which cir- 
cumstances {N) equilibrates. We first note that {N) 
reaches a steady state value only if the interaction energy 
is sufficiently strong. To show this we propagate the ini- 
tial state up to a final time tf = 50 and monitor the 
final value of the expectation of the number of excited 

atoms, (N), and its variance, AA^ = 



{N)\ 



for different interaction strengths V. The correspond- 
ing plot is presented in Fig.[ljb). The thick curves result 
from averaging over the time window A7 < fit < 50 at 
any given value of the interaction strength. Thin lines 
show instantaneous snapshots taken during this time in- 
terval. For V — we find large amplitude oscillations of 
the snapshots which differ strongly from the mean value 
(A'^) « 0.5 L. In this non-interacting regime no relaxation 
into a steady state occurs for the lattice gas subsystem. 
For increasing V, but still with F <C ^2, the instanta- 
neous fluctuations decrease compared to the mean. In 



this regime the dynamics is approximately governed by 
an integrable Hamiltonian as shown in Ref. f23|. Al- 
though (A^) becomes stationary we do not consider this 
to be a thermal state since its fluctuations AA^ are large. 
Moreover, we will later see that a temperature cannot be 
meaningful defined in this weakly interacting regime. 

More interesting for us is the strongly interacting 
regime. It is located in Fig. [TJb) right after the drop 
in AA^ which occurs at V « 2il. In this regime Hamilto- 
nian ll]) is no longer analytically solvable. Here the driv- 
ing term Hq leads to the thermalization of observables 
related to the lattice gas. In particular, we show below 
that the probability p„ to measure n excited particles is 

given by Pn = J^e -P"^ °^ J2e ^^P '^("' ^) ^^P Here 
/3 is an inverse temperature, /3 = (in what follows 

we set the Boltzmann constant to unity, /cs — 1), and 
S{n, e) is an entropy function obtained from the number 
of states with n particles located in a given interval of 
the interaction energy e. 

To prove this we proceed as follows. We define the 
projector P„e = P„ (g) P^, where P„ projects onto the 
subspace containing n particles, and P^ projects onto 
the subspace of states which possesses the interaction 
energy e. As in the numerics we have to deal with a fi- 
nite system, and thus a non-continuous density of states, 
we define a coarse-grained interaction energy: We intro- 
duce energy intervals that are centered around £k = kV 
(with fc = 0, 1, ■■■,L) and have width V. Each of these 
energy intervals contains a large number of (quasi) degen- 
erate states onto which P,. projects. We now transform 
into the interaction picture with the transformation U = 
exp {—i Hj^Q t), where the Hamiltonian for fi = is given 

by H,{t) = WHqU = J2ee' PeHqPe' exp(-lWee't) with 

uj^^r — e' ~ e. The evolution of the density matrix is then 
given by the von-Neumann equation dtp = —i [Hx{t),p\ 
which is equivalent to 



dpjt) 
dt 



= -i [H^{t), p{0)] - fds [H^{t), [H^{s), 
Jo 



We can now to derive a closed equation for the evolution 
of Pne = TrP„£p(i). By generalizing the ideas from Ref. 
we find that this equation is given by 



dp„ 



dt 
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Here we have exploited the fact that due to the strong 
interaction among the particles the couplings between 
subspaces with different particle number n ^ m and 
interaction energy e ^ e' are uncorrelatcd. The mean 
strength of these couplings is contained in the coefficients 
Knm(£,£')- We are now interested in coarse-graining this 
evolution equation over time intervals that are larger 
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FIG. 2: (a) Entropy function for L = 20. Shown is S{n, s) — 
log dim, where dim is the total number of states of the system. 
The interaction energy is coarse grained over the interval Ek = 
[kV- V/2, kV + V/2] . (b) Evolution of (TV) into equilibrium 
for V — 5Q. (c) The logarithm of pe/g{e) is given by the 
dotted line. The red solid line represents a linear fit. 



than il ^. After this procedure one obtains 
dpn 




dt 



(X 



FIG. 3: (a) Inverse temperature as a function of the in- 
teraction strength after a propagation time of nt — 250 
(L = 20). The error bars delimit the 95% confidence in- 
terval, (b) Mean density (A'') /L of the classical lattice gas in 
thermal equilibrium plotted as a function of the interaction 
strength V (red/dash curve). The inverse temperature was 
taken from panel (a). The black solid curve shows the result 
if next-nearest neighbor interactions are accounted for pertur- 
batively. For comparison the mean density which is obtained 
from the quantum calculation is shown (blue circles). 



where the coefficients Umnis) contain the coupling rates 
between subspaces with fixed interaction energy but dif- 
ferent particle number. The quantity dim„(e) is the num- 
ber of states that lie in the interaction energy window 
centered around e and contain n excited particles. Pro- 
vided that Um„(£) 7^ the steady state solution of this 
equation is of the form 

Pne = 6edim„(e) = 6j expS'(n,e). 

The entropy S{n,e) is shown in Fig. [DJa). Since sub- 
spaces with different interaction energy are independent 
and uncoupled the probability to have n excitations is 
given by p„ = J^eP^e = J2e^e expS{n,e). 

What remains is to determine the functional depen- 
dence of the coefficients be on the interaction energy from 
numerical simulations. As an illustration we consider a 
system with L = 20 and V = 5$! whose temporal evolu- 
tion of the mean particle number is shown in Fig. [2][b). 
We study the system a.t t = 250 fi"^ where the steady 
state is well-established. In order to determine the 6e's 
we consider = J^nPne = EnexpS'(n,e) = b^gie) 
where g{e) is the number of states with energy e, and 
study the ratio Pe/g{s). In the Fig. [2}; the logarithm of 
this quantity is plotted as a function of e which is ex- 
cellently fitted by a Boltzmann distribution, and we find 
be cx exp-;9e with {3 = 0.376 f^-^. 

The numerical results strongly suggest that a.tV — 5^1 
the steady state of the system is given by a thermal 
(canonical) distribution of the lattice gas Hamiltonian 
-f^LG with an inverse temperature /3 that is determined 
by the Rabi frequency il. To confirm if this also holds 
for other interaction strengths V we determine the tem- 
perature as a function of V according to the procedure 



outlined above. The corresponding data is shown in Fig. 
[3ja). Below V = 2n we see that the Boltzmann distri- 
bution gives a poor account of the numerical data, i.e., 
the fit of be to a function of the form exp — /3e results in 
huge errors. For higher values of V the fit works, which 
is consistent with the observations of Fig.[ljb). 

We now consider whether the above temperature 
measure gives results that are consistent with the ex- 
pectation values obtained in a statistical mechanics 
way. We do so by calculating the excitation den- 
sity in the thermal ensemble of the lattice gas at 
the estimated temperature and comparing this to the 
actual excitation density from the quantum problem. 
The partition function of the lattice gas is given as 
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thermal average of an observable O is given by (O) 
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In 



S-i(/3,M,^) E{„.}Oexp-/3 
the case of nearest-neighbor interactions on a ring, 
'E.{/3, p,, L) can be calculated straightforwardly via the 
transfer-matrix method 21|, and we use this result ne- 
glecting contributions from interactions over a longer 
range. When necessary these interactions can be in- 
cluded perturbatively, as they are smaller by at least 
a factor of 1/64 than the dominant nearest-neighbor 
contribution. The partition function is given by 
S(/3,/u,L) = + with A± EE i (l-f e/5(^-^)) ± 

In the limit L — > (X), which 
20, the excitation den- 



i^4e/3M + (i_e/3(M"^))^ 
is a good approximation for L 
sity at p, = (i.e., excitation laser on resonance) ob- 
tained from this partition sum takes the simple form 
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FIG. 4: Mean number of excited spins, (A'^) (blue), and 
its variance, (iV^) - (iV)^ (red) for L = 20. At Q.t = 
50 the chemical potential is ramped up linearly from 
to 857. The thick curves result from smoothing the data. 
The inset shows the inverse temperature calculated from the 
fluctuation-dissipation theorem ([2]). The dashed curve is cal- 
culated using the classical lattice gas model at a chemical 
potential ^{t) and the corresponding temperature of the in- 
set. 



In Fig. El^b) we compare the excitation density that 
is obtained form the quantum calculation to the equilib- 
rium value in a classical lattice gas at inverse temperature 
j3 taken from Fig.[3lja). The agreement is good, but sys- 
tematic deviations are evident beyond V ~ A^. These 
deviations are reduced if next-nearest neighbor interac- 
tions are included perturbatively. Other initial states dif- 
ferent from the one considered here are expected to lead 
to a better agreement with the lattice gas prediction, in 
analogy with the observations of Ref. [(| . 

Our findings indicate that the temperature inferred 
from the distribution of the interaction energy is compat- 
ible with a classical thermal ensemble. This temperature 
can be estimated in an alternative way, which corrob- 
orates our claim that the system has reached a thermal 
state with respect to observables of the lattice gas. This is 
done by comparing the response of the system to a change 
of the chemical potential with its spontaneous fluctua- 
tions. In thermal equilibrium the fluctuation-dissipation 
theorem states that the change in the average number of 
excitations is given by its variance, i.e.. 
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the constant of proportionality being the inverse temper- 
ature. Experimentally such fluctuations have been mea- 
sured for example in frozen Rydberg gases where 
they are used to characterize the relative importance of 
We are thus confident that the method 
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interactions 

we outline below will turn out to be a valuable tool for 
the experimental determination of the temperature. 

In Fig. 2] we show a numerical simulation of the envis- 
aged experiment. We start the Hamiltonian evolution of 



the system with = 0, keeping it that way until VLt = 50 
for the system to equilibrate. We then linearly ramp up 
the chemical potential to /imax = Sfl during a time in- 
terval of length 100 r2~^. Experimentally this is done by 
decreasing the detuning of the excitation laser which, as 
expected, leads to an increase of the mean number of 
excited particles. The parameters are chosen such that 
no saturation of the excitation number occurs. We de- 
termine the fluctuations {AN)"^^ = {N^)^^ - (^)?h' and 
obtain the inverse temperature from Eq. ([2|), using that 

d^{N{t))ti, = [^]-'5t(iV(<))th. The (time-dependent) 
temperature is plotted in the inset of Fig. 01 At the begin- 
ning of the ramp (fit = 50) the density does not respond 
immediately to the change of ^. This finite time effect 
causes a spurious initial increase of (3 from zero which 
saturates at fit « 70 at /3 = 0.44 fi"-"^. For consistency, 
this value has to be close to the value predicted from the 
distribution of the interaction energies at the beginning 
of the ramp, which is /? = 0.40 fi~^. Given the strongly 
fluctuating data and the data smoothing the agreement 
is good. The inset of Fig. |4] also shows that with increas- 
ing fit, and hence increasing /i, the inverse temperature 
/3 decreases. The system is heated, showing that the evo- 
lution is not fully adiabatic. For comparison we display 
{N{t))^^ calculated from the classical partition function 
with fi{t) and f3{t) (dashed line). The agreement is ex- 
cellent. 

The fact that both schemes suggested here give consis- 
tent estimates for the temperature is non-trivial. It indi- 
cates that the subsystem of the closed quantum problem 
that corresponds to the lattice gas observables is indeed 
in thermal equilibrium. This shows that aspects of the 
real time dynamics (for long times) of a strongly interact- 
ing and coherently driven quantum many-particle system 
can be understood through classical equilibrium thermo- 
dynamics. 
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